% Approximate Hellinger distance in bivariate example under
% historical-decomposition restriction.

clear variables
close all
clc

rng(121413); % Set random seed

T = 3; % Length of time series
n = 1000; % Number of gridpoints over which to evaluate likelihood
k = 1; % Period in which first structural shock is positive
K = 1e4; % Monte Carlo sample size for approximating Hellinger distance
Keps = 1e6; % Monte Carlo sample size for estimating probabilities

% Set true values of structural parameters.
A0 = [1 0.5; 0.2 1.2];
% Matrix of impulse responses.
B = inv(A0);

% Compute reduced-form parameters
% Lower-triangular Cholesky factorisation of innovation covariance matrix.
Sigmatr = chol(inv(A0)*inv(A0)','lower'); 
Sigma = Sigmatr*Sigmatr'; % Covariance matrix of VAR innovations.
sig11 = Sigmatr(1,1);
sig21 = Sigmatr(2,1);
sig22 = Sigmatr(2,2);

% Back out corresponding values of theta in orthogonal reduced form.
thetaTrue = asin(A0(1,2)*sig22);

% Construct grid for theta over which to evaluate Hellinger distance.
thetaGrid = linspace(-pi+eps,pi-eps,n)';

%% Compute Hellinger distance for unconditional likelihood.

hellinger = zeros(n,1);

for kk = 1:K
    
    % Draw a time series of structural shocks and check whether narrative
    % restrictions hold (the first structural shock in the kth period is
    % positive and the most important contributor to the observed change in
    % the first variable).
    epst = randn([T,2]);
    flagHD = (abs(B(1,1)*epst(k,1)) >= abs(B(1,2)*epst(k,2)));
    flagSR = (epst(k,1) >= 0);

    % Use shocks and structural parameters to generate data.
    yt = (A0\epst')';
   
    unLik = zeros(n,1);
    
    for ii = 1:n % For each value of theta

        theta = thetaGrid(ii);

        % Value of first structural shock given value of theta and 
        % realisation of data in period k.
        eps1k = (1/(sig11*sig22))*(sig22*yt(k,1)*cos(theta) ...
            + (sig11*yt(k,2)-sig21*yt(k,1))*sin(theta));

        % Check whether narrative restriction and sign normalisations are 
        % satisfied at value of theta.

        if cos(theta) >= 0 && sig22*cos(theta) - sig21*sin(theta) >= 0 ...
                && (eps1k*flagSR -eps1k*(1-flagSR) >= 0)

            % Value of second structural shock given value of theta and 
            % realisation of data in period k.        
            eps2k = (1/(sig11*sig22))*((sig11*yt(k,2)-sig21*yt(k,1))*cos(theta) ...
                - sig22*yt(k,1)*sin(theta));

            if (flagHD == 1 && abs(sig11*cos(theta)*eps1k) >= ...
                    abs(-sig11*sin(theta)*eps2k)) || ...
               (flagHD == 0 && abs(sig11*cos(theta)*eps1k) < ...
                    abs(-sig11*sin(theta)*eps2k))
                                     
                % Add contribution to Hellinger distance.
                hellinger(ii) = hellinger(ii) + 1/K;

            end

        elseif cos(theta) < 0 && sig22*cos(theta) - sig21*sin(theta) >= 0 ...
                && (eps1k*flagSR -eps1k*(1-flagSR) >= 0)

            % Value of second structural shock given value of theta and 
            % realisation of data in period k. 
            eps2k = (1/(sig11*sig22))*(sig22*yt(k,1)*sin(theta) ...
                -(sig21*yt(k,1) + sig11*yt(k,2))*cos(theta));

            if (flagHD == 1 && abs(sig11*cos(theta)*eps1k) >= ...
                    abs(sig11*sin(theta)*eps2k)) || ...
               (flagHD == 0 && abs(sig11*cos(theta)*eps1k) < ...
                    abs(sig11*sin(theta)*eps2k))
                 
                % Add contribution to Hellinger distance.
                hellinger(ii) = hellinger(ii) + 1/K;

            end

        end

    end   
    
end

HDsq = 2*(1 - hellinger); % Squared Hellinger distance

%% Compute Hellinger distance for conditional likelihood.

% Draw structural shocks from bivariate standard normal distribution.
epsk = randn([Keps,2]);

% Compute probability that structural shocks satisfy narrative restriction 
% at each value of theta on grid.

PrEps = zeros(n,1);

for ii = 1:n

    PrEps(ii) = compPrEps(Sigmatr,thetaGrid(ii),epsk);

end

% Compute Hellinger distance at each value of theta on grid.

hellingerC = zeros(n,1);

for kk = 1:K
    
    % Draw a time series of structural shocks such that the first 
    % structural shock in the kth period is the most important contributor 
    % to the observed change in the first variable.
    epst = randn([T,2]);
    flag = 0;

    while flag == 0

        epst(k,1) = randn;
        flag = ((abs(B(1,1)*epst(k,1)) >= abs(B(1,2)*epst(k,2))) & ...
            epst(k,1) >= 0);

    end

    % Use shocks and structural parameters to generate data.
    yt = (A0\epst')';
    
    Lik = zeros(n,1);
    
    for ii = 1:n % For each value of theta

        theta = thetaGrid(ii);

        % Value of first structural shock given value of theta and 
        % realisation of data in period k.
        eps1k = (1/(sig11*sig22))*(sig22*yt(k,1)*cos(theta) ...
            + (sig11*yt(k,2)-sig21*yt(k,1))*sin(theta));

        % Check whether narrative restriction and sign normalisations are 
        % satisfied at value of theta.

        if cos(theta) >= 0 && sig22*cos(theta) - sig21*sin(theta) >= 0 ...
                && eps1k >= 0

            % Value of first structural shock given value of theta and 
            % realisation of data in period k.        
            eps2k = (1/(sig11*sig22))*((sig11*yt(k,2)-sig21*yt(k,1))*cos(theta) ...
                - sig22*yt(k,1)*sin(theta));

            if abs(sig11*cos(theta)*eps1k) >= abs(-sig11*sin(theta)*eps2k)

                % Compute conditional likelihood.
                Lik(ii) = exp(- T*log(2*pi) - (T/2)*log(det(Sigma)) ...
                    - 0.5*trace((yt/Sigma)*yt') - log(PrEps(ii)));

            end

        elseif cos(theta) < 0 && sig22*cos(theta) - sig21*sin(theta) >= 0 ...
                && eps1k >= 0

            % Value of first structural shock given value of theta and 
            % realisation of data in period k. 
            eps2k = (1/(sig11*sig22))*(sig22*yt(k,1)*sin(theta) ...
                -(sig21*yt(k,1) + sig11*yt(k,2))*cos(theta));

            if abs(sig11*cos(theta)*eps1k) >= abs(sig11*sin(theta)*eps2k)

                 % Compute conditional likelihood.
                Lik(ii) = exp(- T*log(2*pi) - (T/2)*log(det(Sigma)) ...
                    - 0.5*trace((yt/Sigma)*yt') - log(PrEps(ii)));                                

            end

        end

    end
    
    % Compute likelihood at true value of theta.
    PrEpsTrue = compPrEps(Sigmatr,thetaTrue,epsk);
    LikTrue = exp(- T*log(2*pi) - (T/2)*log(det(Sigma)) ...
                    - 0.5*trace((yt/Sigma)*yt') - log(PrEpsTrue));  

    % Compute Hellinger distance.
    hellingerC = hellingerC + sqrt(Lik./LikTrue)./K;
    
end

HDCsq = 2*(1 - hellingerC); % Squared Hellinger distance

% Plot Hellinger distances.
figure;
h1 = plot(thetaGrid,HDsq,'color',[255 0 0]./255,'LineWidth',2);
hold on;
h2 = plot(thetaGrid,HDCsq,'color',[255 0 255]./255,'LineWidth',2);
h3 = line([thetaTrue thetaTrue],[0 2],'color',[76 0 153]./255,'LineWidth',2);
xlim([-pi pi]);
xticks([-pi -pi/2 0 pi/2 pi]);
xticklabels({'-\pi','-\pi/2','0','\pi/2','\pi'});
xlabel('\theta');
ylim([0 2]);
title('Historical Decomposition');
legend([h1 h2 h3],{'Unconditional','Conditional','$$\theta_{0}$$'},...
    'interpreter','LaTeX','location','SouthWest','FontSize',14);
legend boxoff;
Ax = gca;
Ax.FontSize = 14;
cd('Figures');
print('HistoricalDecompositionHellinger','-depsc');
print('HistoricalDecompositionHellinger','-dpng');
cd ..

cd('Results');
save('HistoricalDecompositionHellinger_results.mat');
cd ..

%%
function prob = compPrEps(Sigmatr,theta,epsk)
% Compute probability that shocks satisfy narrative restrictions.
% Inputs:
% - Sigmatr: lower-triangular Cholesky factor of Sigma.
% - theta: structural parameter
% - epsk: matrix containing Monte Carlo sample of structural shocks

sig11 = Sigmatr(1,1);
sig21 = Sigmatr(2,1);
sig22 = Sigmatr(2,2);

if cos(theta) >= 0 && sig22*cos(theta) - sig21*sin(theta) >= 0
    
    prob = mean((abs(sig11*cos(theta)*epsk(:,1)) >= ...
        abs(-sig11*sin(theta)*epsk(:,2))) & epsk(:,1) >= 0);
    
elseif cos(theta) < 0 && sig22*cos(theta) - sig21*sin(theta) >= 0
    
    prob = mean((abs(sig11*cos(theta)*epsk(:,1)) >= ...
        abs(sig11*sin(theta)*epsk(:,2))) & epsk(:,1) >= 0);
    
else 
    
    prob = 0;
    
end

end